
#######################################################################################
# File: MESMLMLogit-AR0.R                                                             
# Author: Xun PANG                                                             
# Creation Date: 4/5/2014                                                            
# Description: 	MCMC simulation for estimating binary Logit MESMLM-AR(0) model
#               based on the research paper, "Varying Responses to Common Shocks and
#               Complex Cross-Sectional Dependence: Dynamic Multilevel Modeling with
#               Multifactor Error Structures for Time-Series Cross-Sectional Data "
#
# Usage:         MulFacLogit.AR0 <- function(y,X1, Wi, Ai, Ut, TP, Unit.Index, ...)
#    
# Arguments:
#         y     response variable, dichotomous (0/1)
#        X1     matrix of covariates with fixed effects
#        Wi     matrix of covariates with subject-varying effects
#        St     matrix of covariates with time-varying effects
#        Ai     matrix of covariates explaining the subject-varying effects
#        TP     number of unobserved factors
# Unit.Index    index of subjects
# Time.Index    index of time periods
# timeint.add   add an time-specific intercept? "TRUE" or "FALSE". Default is "FALSE"
# unitint.add   add an unit-specific intercept? "TRUE" or "FALSE". Default is "FALSE"
#         m     number of iterations to be returned
#    burnin     number of initial iterations which are discarded
#    priors     beta0--- mean vector of the multivariate normal distribution of beta,
#                        all the fixed effects (beta in the reduced form of the model)
#               B0--- variance matrix of the multivariate normal distribution of beta,
#                        all the fixed effects (beta in the reduced form of the model)
#               d0, D0 --- degree of freedom and scale matrix of the Wishart prior
#                          on b_i, the subject-level residual
#               e0, E0 --- degree of freedom and scale matrix of the Wishart prior
#                          on c_t, the time-level residual
# initial values
#               init.phi --- vector of initial values of the autoregressive coef.
#               init.beta, init.b, and init.c --- initial values of the fixed and
#                   random coefficients
#               init.ystar --- initial values of the latent response variable
# tracking      every "tacking" iterations, return the information about how many
#               iterations in total have been done
#  monitor      a string contains the names of parameters whose MCMC output are
#               will be returned. The string has to be a subset of
#               ("rho", "beta", "bi", "ct", "D", "E", "ystar", "u") which is the default
#     thin      thin the chain by recording each "thin" iterations
#                                       
###########################################################################################################
# Call functions
source("hidden.function.R")
## Define SeveralFunctions for Rejection Sampling


  ################################################################
  # function to sample from truncated logistic distribution
  ################################################################

  # x is the observed dichotomous variable
  # mu: mean of logistic distrbution
  # sigma: scale parameter
  BTlogisticlog <- function(x, mu, sigma=1){
    if(x==0) {
       u <- runif(1)
       xstar <- qlogis(log(u)+plogis(0, mu, sigma, log.p=T), mu, sigma, log.p=T)
    } else{
       u <- runif(1)
       xstar <- -qlogis(log(u)+plogis(0, -mu, sigma, log.p=T), -mu, sigma, log.p=T)
    }
    xstar
  }
  ################################################################
  # function to sample from truncated normal distribution
  ################################################################

  # x is the observed dichotomous variable
  # mu: mean of normal distrbution
  # sigma: standard deviation
  BTnormlog <- function(x, mu, sigma=1){
    if(x==0) {
       u <- runif(1)
       xstar <- qnorm(log(u)+pnorm(0, mu, sigma, log=T), mu, sigma, log=T)
    } else{
       u <- runif(1)
       xstar <- -qnorm(log(u)+pnorm(0, -mu, sigma, log=T), -mu, sigma, log=T)
    }
    xstar
  }
   # Two functions for the rejecting sampling in two partitionalizations of the sample space
   rightmost.interval <- function(U, lambda, mm){
    ZD <- 1
    XD <- exp(-0.5*lambda)
    j <- 0
    for (d in 1:mm){
      j <- j+1
      ZD <- ZD-(j+1)^2*XD^((j+1)^2-1)
      if (ZD>U) return(TRUE)
      j <- j+1
      ZD <- ZD+(j+1)^2*XD^((j+1)^2-1)
      if (ZD<U) return(FALSE)
    }
    return(FALSE)
  }

  leftmost.interval <- function(U, lambda, mm){
    HD <- 0.5*log(2)+2.5*log(pi)-2.5*log(lambda)-pi^2/(2*lambda)+0.5*lambda
    LU <- log(U)
    ZD <- 1
    XD <- exp(-pi^2/(2*lambda))
    KD <- lambda/pi^2
    j <- 0
    for (d in 1:mm){
      j <- j+1
      ZD <- ZD-KD*XD^(j^2-1)
      WD <- HD+log(ZD)
      if (WD>LU) return(TRUE)
      j <- j+1
      ZD <- ZD+(j+1)^2*XD^((j+1)^2-1)
      WD <- HD+log(ZD) 
      if (WD<LU) return(FALSE)
    }
    return(FALSE)
  }

  ## Function for sampling from generalized inverse gaussian distribution
  ## n is the number of draws, and jj is one of the parameters
  dGIG <- function(n, jj){
       draw <- c()
      for (i in 1:n){
        R <- (rnorm(1))^2
        RR <- 1+(R-sqrt(R*(4*jj+R)))/(2*jj)
        RRR <- 1/(1+RR)
        u <- runif(1)
        if (u < RRR){
          draw[i] <- jj/RR
        }else{
          draw[i] <- jj*RR
        }
      }
        draw
     }

  ###################################
  # Function for MCMC updating
  ###################################
  MulFacLogit.AR0 <- function(y,X1, Wi, Ai, Ut, TP, Unit.Index, Time.Index,
                              m=10000, burnin=5000,  init.ystar=0, init.beta=0,
                              unitint.add=0,timeint.add=0, init.b=0, init.c=0, init.theta=0,
                              init.gamma=0, beta0, B0, D0, d0, E0, e0,  mm=100, 
                              monitor=c( "beta", "bi", "ct", "D", "E", "ystar", "theta", "gamma"),
                              thin=1, tracking){
    
    #****************** Arrange Data for MCMC updating  **********************#
     
     # To use the RegressionData for this model specification, the arguments of
     #  Ft and St are set to be certain values since the data inpute does not contain
     #  information for Ft and St
     Ft <- "NULL"
     St <- as.matrix(rep(1, length(y)))
     # Call RegressionData() and rearrange the data to fit the reduced model
     NewData <- RegressionData(yob=y, Unit=Unit.Index, Time=Time.Index, Fixed=X1,
                               UnitRandom=Wi, TimeRandom=St,
                               UnitPred=Ai, TimePred=Ut)

    Y <- NewData[[3]]
    entry1 <- 1: length(Y)
    newid1 <- NewData[[1]]
    newid2 <- NewData[[2]]
    id11 <- Index(newid1)     #the unit diex
    uniqunit <- as.numeric(table(id11))   # how many time periods for each unit
    jj <- which(uniqunit<= nlag)
    if (length(jj)>0){						
      entry2 <- entry1[id11==jj]
      Y <- Y[-entry2]
      X <- as.matrix(NewData[[4]][-entry2,])
      WP <- as.matrix(NewData[[5]][-entry2,])
    }else{
      X <- as.matrix(NewData[[4]])
      WP <- as.matrix(NewData[[5]])
    }
    WPP <- ncol(WP)
    KW <- ncol(Ut)
    FB <- ncol(X)
    W <- as.matrix(cbind(WP, X[, (FB-KW+1):FB]))
    if (unitint.add==TRUE){
      W <- cbind(1, W)
    }
   N <- length(Y)
     
   # Check the dimension of the parameter vectors
     
   #beta: fixed parameter vector 
   PFB <- length(beta0)
   if (FB!=PFB){
      cat("Length of the prior on fixed effects is incorrect.\n")
      cat("Length of fixed effect should be ", FB, "\n", sep=" ")
      stop("Respecify beta0 and B0, call GLMMARp.Binary() again. \n") 
    }
   # bi: random fixed parameter vector 
   RB <- ncol(W)
   if(is.matrix(D0)){ 
     PRB <- ncol(D0)-TP
   }else{
     PRB <- length(D0)-TP
   }
   if (RB!=PRB){
       cat("Dimension of the prior on unit random effects is incorrect.\n")
       cat("Length of D0 should be ", RB+TP, "*", RB+TP,  "\n", sep=" ")
       stop("Respecify D0, call GLMMARp.Binary() again. \n") 
    }

   # in case of unbalanced TSCS data structure,
   # create new indices to accomendate such a situation 						
   entry <- seq(1: N)                   # the global length of observations
   if (length(jj)>0){
     # unit id repeating t times t=# observations
     id1 <- Index(id11[-entry2])                
     # year id repeating n times n=# unit in each time period
     id2 <- (newid2[-entry2])-min(newid2[-entry2])+1        
     GU <- max(id1)                       # how many units
     GT <- max(id2)                       # how many time periods in total
     id3 <- id1
     GU <- length(unique(id3))
     n=1
     while (n <= length(unique(id3))){
       rows <- entry[which(id3==n)]                   
       nT <- length(rows)
       TT <- id2[rows]
       if (nT!=1){
         for (j in 2:nT){
          GHK <- TT[j]-TT[j-1]
          if (GHK >1){
            PL <- rows[j]
            id3[PL:N] <- id3[PL:N]+1
          }
        }
       }
       n <- n+1
     }

    entry3 <- 1: length(Y)
    uniqunit <- as.numeric(table(id3))
    jj <- which(uniqunit<= nlag)
    entry4 <- entry1[id3==jj]
    if(length(entry4)>0){
      Y <- Y[-entry4]
      X <- as.matrix(X[-entry4,])
      W <- as.matrix(W[-entry4,])
      id1 <- Index(id1[-entry4])
      id2 <- Index(id2[-entry4])
      id3 <- Index(id3[-entry4])
    }
   }else{
     id1 <- id11
     id2 <- Index(newid2)
   }
    id1 <- Index(id1)                 # unit id repeating t times t=# observations
    id2 <- (id2)-min(id2)+1        # year id repeating n times n=# unit in each time period
    id3 <- id1						
    GU <- max(id1)                       # how many units
    GT <- max(id2)
    GUU <- length(unique(id3))


     #############################
     # Storage and Initial Values
     #############################

    total.return <- m+burnin

   if ("beta"%in% monitor){ 
     beta <- matrix(NA, nrow=total.return, ncol=FB)
     beta[1,] <- init.beta
   }else{
     beta <- matrix(init.beta, nrow=1, ncol=FB)
   }
   
   if ("bi"%in% monitor){ 
     bbb <- matrix(NA, nrow=total.return, ncol=RB*GU)
     bbb[1,] <- matrix(init.b, nrow=1, ncol=RB*GU)
   }else{
     bbb <- matrix(init.b, nrow=1, ncol=RB*GU)
   }

   if (timeint.add==1) {
     TPP <- TP+1
   }else{
     TPP <- TP
   }
   if ("ct"%in% monitor){ 
    ccc <- matrix(NA, nrow=total.return, ncol=GT*TPP)
    ccc[1,] <- matrix(init.c, nrow=1, ncol=GT*TPP)
   }else{
    ccc <- matrix(init.c, nrow=1, ncol=GT*TPP)
   }
   EDim <- (TPP^2-TPP)/2+TPP
   if ("E"%in% monitor){
     EEE <- matrix(NA, nrow=total.return, ncol=EDim)
     EEE[1,] <- matrix(0, nrow=1, ncol=EDim)
   }else{
     EEE <- matrix(0, nrow=1, ncol=EDim)
   }
   
   DDim <- ((RB+TP)^2-(RB+TP))/2+RB+TP
   if ("D"%in% monitor){
     DDD <- matrix(NA, nrow=total.return, ncol=DDim)
     DDD[1,] <- matrix(0, nrow=1, ncol=DDim)
   }else{
     DDD <- matrix(0, nrow=1, ncol=DDim)
   }

  if ("theta"%in% monitor){
     theta <- matrix(NA, nrow=total.return, ncol=TPP*GT)
     theta[1,] <- matrix(init.theta, nrow=1, ncol=TPP*GT)
   }else{
     theta <- matrix(init.theta, nrow=1, ncol=TPP*GT)
   }

   if ("gamma"%in% monitor){
     gamma <- matrix(NA, nrow=total.return, ncol=TP*GU)
     gamma[1,] <- matrix(init.gamma, nrow=1, ncol=TP*GU)
   }else{
     gamma <- matrix(init.gamma, nrow=1, ncol=TP*GU)
   }
   
   if ("ystar"%in% monitor){ 
     ystarm <- matrix(NA, nrow=total.return, ncol=N)
     ystarm[1,] <- matrix(0, nrow=1, ncol=N)
   }else{
     ystarm <- matrix(0, nrow=1, ncol=N)
   }
							  
 ##################################################
  # Calculate those quantities which will
  # not be updated in MCMC iterations
  ################################################## 
     invB <- pseudoinverse(B0)
     PBbeta <- invB%*%beta0
     En <- (e0+GT)/2
     En2 <- En*2
     if(is.matrix(E0)){
       invE <- pseudoinverse(E0)
     }
     Dn <- (e0+GT)/2
     Dn2 <- Dn*2
     if(is.matrix(D0)){
       invD <- pseudoinverse(D0)
     }
     cat("Data are arranged. MCMC begins.\n")
     
  ################################################
  ## MCMC updating loop begins here
  ################################################

  #Tracking the process of MCMC simulation
    mcmc <- (m+burnin)*thin
    for(g in 2:(mcmc)) {
      if (tracking>0){
       if ((g%% tracking ==0)| g==(m*thin)){
          cat("g=",g,"\n")
        }
     }
             
   # The current values of those parameters

   if (g==2){
     current.bi <- matrix(bbb[1,], ncol=RB)      # current values of ci and beta3i
     current.theta <- matrix(theta[1,], ncol=TPP)      # current values of thetat
     current.gamma <- matrix(gamma[1,], ncol=TP)      # current values of gammai
     current.beta <- beta[1,]                         # current values of beta
     current.ystar <- ystarm[1,]
     current.lambda <- rep(1, N)                      # the auxilary variable
   }else{
     current.bi <- new.bi                            # current values of ci and beta3i
     currect.theta <- new.theta                      # current values of theta.t
     current.gamma <- new.gamma                      # current values of gamma.i
     current.beta <- new.beta                        # current values of beta
     current.ystar <- new.ystar                      # the current value of ystar
     current.lambda <- new.lambda
   }
  
  ###################################################
  # Block 1:  Covariance Matrices (random errors)
  ###################################################  

   if (TPP==1){ # if S only contains 1 element (most likely a constant)
     Ei<-rgamma(1,En, rate=(E0+sum(current.theta^2))/2)
     Ei <- as.matrix(Ei)
   }else{
     OEin <- invE+sumtrans(Mat=current.theta)
     sE<-rWishart (1, En2, pseudoinverse(OEin))
     Ei <- Wishart(element=sE, col=TPP)
   }

       
   current.unit <- cbind(current.bi, current.gamma)

   if (ncol(current.unit)==1){ # if W only contains 1 element (most likely a constant)
     Di<-rgamma(1,Dn, rate=(D0+sum(current.unit^2))/2)
     Di <- as.matrix(Di)
   }else{
     ODin <- invD+sumtrans(Mat=current.unit)
     sD<-rWishart (1, Dn2, pseudoinverse(ODin))
     Di <- Wishart(element=sD, col=(RB+TP))
  }

    ###################################################
    # Block 2:  data augmentation for Y (and prepare for beta )
    ###################################################
       
    betaterm1 <- matrix(0, nrow=FB, ncol=FB)     # for posterior covariance of beta
    betaterm2 <- rep(0, FB)                      # for posterior mean of beta 
    new.ystar <- rep(NA, N)
    Colum <- (RB+TP)
    new.unit2 <- matrix(NA, nrow=GU, ncol=Colum)
       
    # Big loop for data augmentation and values for beta 
    for (n in 1: GU){                                 # loop for each unit
      rows <- entry[which(id1==n)]                    # the obserbations of unit n
      lambdan <- current.lambda[rows]
      nT <- length(rows)                              # length of those observations
      yi <- Y[rows]                                   # Observed responses of n
      ystari <- current.ystar[rows]                   # current latent responses of n
      Xi <- as.matrix(X[rows,])
      Wi <- as.matrix(W[rows,])
      Si <- current.gamma[n,]
      if (timeint.add==1){
        Si <- c(1, Si)
      }
      bbbn <- current.bi[n,]
      TT <- id2[rows]
      LT <- max(TT)-min(TT)+1
      KIU <- min(TT):max(TT)
      AS <- KIU%in% TT

      Tccin <- as.matrix(current.theta[TT,])
      rant <- c()
      for (l in 1:nT){
        rant[l] <- Si%*%as.matrix(Tccin[l,])
      }

    # Updata ystar
    mean.ystar <- Xi%*% cbind(current.beta)+Wi%*%cbind(bbbn)+rant
    ystar <- sapply(c(1:nT),
             function(r){BTnormlog(x=yi[r], mu=mean.ystar[r], sigma=sqrt(lambdan[r]))})
 
    new.ystar[rows] <- ystar

   #*************************************************************#
    #  Begin to calculate the quantities for updating {ct,f} and gamma
   inv.lambdan <-pseudoinverse(diag(lambdan)) 
   if (n==1){
      if (timeint.add==1) Tccin <- as.matrix(Tccin[,-1])
      NW <- cbind(Wi, Tccin)
      Cov.b <- pseudoinverse(Di + t(NW)%*%inv.lambdan%*%NW)
      mean.b <- Cov.b%*% t(NW)%*%inv.lambdan%*%(ystar-Xi%*% cbind(current.beta))
      u1 <- mean.b[1:RB]
      sigma11 <- Cov.b[1:RB, 1:RB]
      x2 <- rep(1, TP)
      if (length(u1)==1){
        bb.draw <- rnorm(1, u1, sqrt(sigma11))
        b.drawn <- as.matrix(c(bb.draw, x2))
      }else{
        bb.drawn <- mvrnorm(1, mu=u1, Sigma=sigma11)
        b.drawn <- as.matrix(c(bb.drawn, x2))
      } 
    }else{
      if (timeint.add==1) Tccin <- as.matrix(Tccin[,-1])
      NW <- cbind(Wi, Tccin)
      Cov.b <- pseudoinverse(Di + t(NW)%*%inv.lambdan%*%NW)
      mean.b <- Cov.b%*% t(NW)%*%inv.lambdan%*%(ystar-Xi%*% cbind(current.beta))
      if (length(mean.b)==1){
        b.drawn <- as.matrix(rnorm(1, mean.b, sqrt(Cov.b)))
      }else{
        b.drawn <- mvrnorm(1, mu=mean.b, Sigma=Cov.b)
     }
   }
   new.unit2[n,] <- as.vector(b.drawn)
      
   nbbbn <- as.matrix(b.drawn[-c((Colum-TP+1):Colum)])
   ngamma <- as.matrix(b.drawn[c((Colum-TP+1):Colum)])
   nrant <- c()
   for (l in 1:nT){                                # Using loop to do it
     nrant[l] <- Tccin[l,]%*%ngamma
   }
           
   # calculate BB1 and BB2 for beta
    
    betaterm1 <- betaterm1+t(Xi)%*%inv.lambdan %*%Xi
    betaterm2 <- betaterm2+t(Xi)%*%inv.lambdan%*% (ystar-Wi%*% cbind(nbbbn)-nrant)
  }

       
     ###################################
     ##Block 3: Beta and {b_i, \gamma_i}
     ###################################
      # 1) {b_i, \gamma_i}
     new.bi <- as.matrix(new.unit2[, -c((Colum-TP+1):Colum)])
     new.gamma <- as.matrix(new.unit2[,  c((Colum-TP+1):Colum)])

     # 2) \beta
     B1<-pseudoinverse(betaterm1+invB)
     beta1<-t(B1%*%(betaterm2+PBbeta))
     new.beta<-mvrnorm(1, beta1, B1)
        
    #######################################
    # Block 4: Time-Specific Random Effect
    ######################################
      
  new.theta <-matrix(NA, nrow=GT, ncol=TPP)
  for (j in 1: GT){
    row1 <- entry[which(id2==j)]    # Those entries with t=j
    lambdat <- current.lambda[row1]
    inv.lambdat <-pseudoinverse(diag(lambdat))
    H <- length(row1)               #  #observations at t period
    iden1 <- id1[row1]              # The units at t period
    tttn <- as.matrix(new.bi[iden1,])         # The random unit-specific effects of
    yt1 <- new.ystar[row1]
    Xt1 <- as.matrix(X[row1,])
    Wt1<- as.matrix(W[row1,])
    St1 <- as.matrix(new.gamma[iden1,])
    if(timeint.add==1){
      St1 <- cbind(1, St1)
    }
            
    term7 <- c()
    for (q in 1:H){
      term7[q] <- Wt1[q,]%*% cbind(tttn[q,])
    }
    cov.S <- pseudoinverse(Ei+t(St1)%*%inv.lambdat%*%St1)
    mean.S <- cov.S%*% t(St1)%*%inv.lambdat%*%(yt1-term7-Xt1%*%cbind(new.beta))
    if (length(mean.S)==1){
       supd <- as.matrix(rnorm(1, mean.S, sqrt(cov.S)))
    }else{
       supd <- rmvnorm(1, mean.S, cov.S)
    }
    new.theta[j,] <- as.vector(supd)
  }
 
   #########################################################
   ## Block 5: update lambda using rejecting sampling
   #########################################################

    new.lambda<- c()
    for (n in 1: GU){                                 # loop for each unit
      rows <- entry[which(id1==n)]                    # the obserbations of unit n
      nT <- length(rows)                              # length of those observations
      ystari <- new.ystar[rows]                   # current latent responses of n
      Xi <- as.matrix(X[rows,])
      Wi <- as.matrix(W[rows,])
      Si <- new.gamma[n,]
      if (timeint.add==1){
        Si <- c(1, Si)
      }
      bbbn <- new.bi[n,]
      TT <- id2[rows]
      LT <- max(TT)-min(TT)+1
      KIU <- min(TT):max(TT)
      AS <- KIU%in% TT

      Tccin <- as.matrix(new.theta[TT,])
      rant <- c()
      for (l in 1:nT){
        rant[l] <- Si%*%as.matrix(Tccin[l,])
      }

    rsquare <- (ystari-Xi%*% cbind(new.beta)+Wi%*%cbind(bbbn)+rant)^2
     lambdan.new <- c()
     for (i in 1:nT){
       OK <- FALSE
       while (!OK){
        # draw a proposal
        prop.lambda <- c()
        jj <- sqrt(rsquare[i])
        R <- (rnorm(1))^2
        RR <- 1+(R-sqrt(R*(4*jj+R)))/(2*jj)
        RRR <- 1/(1+RR)
        u <- runif(1)
        if (u < RRR){
          prop.lambda <- jj/RR
        }else{
          prop.lambda <- jj*RR
        }
        
        # accept/reject
        U <- runif(1)
          if (prop.lambda > 4/3){
             OK <- rightmost.interval(U, prop.lambda, mm)
          }else{
             OK <- leftmost.interval(U, prop.lambda, mm)
          }
        }
        lambdan.new[i] <- prop.lambda
     }
      new.lambda[rows] <- lambdan.new
    }
       
   ############################################
   #Store the updated values
   ############################################     
       
  if ((g%%thin ==0)| g==m){
         gg <- g%/%thin
         if ("ystar"%in% monitor){
           ystarm[gg,] <- new.ystar
         }

       
         if ("bi" %in% monitor){
           bbb[gg,] <- as.vector(new.bi)
         }
        
         
         if ("gamma" %in% monitor){
           gamma[gg,] <- as.vector(new.gamma)
         }
        
         if ("theta"%in% monitor){
           theta2 <- as.vector(new.theta)
           theta[gg,] <- theta2
         }
         
         if ("beta"%in% monitor){
           beta[gg,] <- new.beta
         }
       
         if ("D"%in% monitor){
           if ((RB+TP)==1){
             DDD[gg,] <-  as.vector(Di)
           }else{
             DDD[gg,] <-  sD
           }
          }
       
         if ("E"%in% monitor){
           if (TPP==1){
           EEE[gg,] <- as.vector(Ei)
         }else{
           EEE[gg,] <- sE
           }
         }
       }
    }

     ###############################
     # Return the MCMC results
     ###############################  
  if ("beta"%in% monitor){ 
     posterior.beta <- beta[-c(1:burnin),]
  }else{
     posterior.beta <- beta
   }
 if ("bi"%in% monitor){ 
     posterior.unit.res <- bbb[-c(1:burnin),]
  }else{
     posterior.unit.res <- bbb
  }
 if ("gamma"%in% monitor){ 
     posterior.timeunit.res <- gamma[-c(1:burnin),]
  }else{
     posterior.timeunit.res <- gamma
  }
  if ("theta"%in% monitor){ 
     posterior.timetime.res <- theta[-c(1:burnin),]
  }else{
     posterior.timetime.res <- theta
  } 
  if ("E"%in% monitor){ 
     posterior.time.cov <- EEE[-c(1:burnin),]
  }else{
     posterior.time.cov <- EEE
  }    
  if ("D"%in% monitor){ 
     posterior.unit.cov <- DDD[-c(1:burnin),]
  }else{
     posterior.unit.cov <- DDD
  }
   if ("ystar"%in% monitor){ 
     posterior.data <- ystarm[-c(1:(burnin/thin)),]
  }else{
     posterior.data <- ystarm
  }

    return(list(Beta=posterior.beta,  bi=posterior.unit.res,
                gamma=posterior.timeunit.res, theta=posterior.timetime.res,
                var.theta=posterior.time.cov, Cov.betai=posterior.unit.cov,
                data.augmented=posterior.data))
   }
